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We consider magnetic flux penetration in a superconducting film with a concave corner. Unlike 
convex corners, where the current flow pattern is easily constructed from Bean's critical state model, 
the current flow pattern at a concave corner is highly nontrivial. To address the problem, we do 
a numerical flux creep simulation, where particular attention is paid to efficient handling of the 
non-local electrodynamics, characteristic of superconducting films in the transverse geometry. We 
find that the current stream lines at the concave corner are close to circular, but the small deviation 
from exact circles ensure that the electric field is finite and continuous. Yet, the electric field is, 
as expected, very high at the concave corner. At low fields, the critical state penetration is deeper 
from the concave corner than from the straight edges, which is a consequence of the electrodynamic 
non-locality. A magneto-optical experiment on a YBa2CusO x displays an almost perfect match 
with the magnetic flux distribution from the simulation, hence verifying the necessity of including 
electrodynamic non-locality in the modelling of superconducting thin films. 

PACS numbers: 74.25.Ha, 74.25.Op 



I. INTRODUCTION 

The macroscopic magnetic properties of type-II super- 
conductors can to a large extent be described by the crit- 
ical state model, first formulated by Bean [TJ [2] . It is 
based on the assumption that the magnitude of the cur- 
rent cannot exceed the critical current density j c , thus 
limiting the ability of the superconductor to carry trans- 
port current or to shield applied magnetic fields. The 
original critical state model was formulated for bulk sam- 
ples, but it has later been extended to thin films, where 
additional shielding currents j < j c flow in the region be- 
yond the flux penetration front. Analytical results have 
been found for an infinite strip [3] [4] and circular disk 
[5j [6], while for less symmetric shapes or disconnected 
geometries the critical state has only been determined 
numerically [7]. 

Dynamics is disregarded in Bean's critical state model, 
in the sense that the response to an applied magnetic 
field is instantaneous. In particular, the model does not 
take into account flux creep, a process acting in any real 
type-II superconductor, and is particularly pronounced 
in high-T c superconductors. To model flux creep, time 
must be taken explicitly into consideration. Convention- 
ally, this is done by assuming a highly nonlinear E — J re- 
lation before solving Maxwell's equations. For thin films, 
the main obstacle for an efficient implementation of a nu- 
merical scheme, is the handling of the boundary condi- 
tions, given the nonlocality of the governing equations. 
Brandt has derived solutions for selected geometries, such 
as rectangles [8 , disks and rings [9], or arbitrary shape 
[TO] , based on a matrix inversion method. Although be- 
ing accurate, these solutions scale poorly with system 
size and need 0(N 2 ) operation for each time step, where 
N is the number of discrete points in the spatial grid. 



Better scaling properties can be obtained by the conju- 
gate gradient method, which scales like 0(N 1A ) pT| [T2] 
or a hybrid real space - Fourier space method scaling as 

0(Nlog(N)) Hani. 

One remarkable property of the Bean model is that 
the current stream line pattern can be drawn simply by 
adding lines with constant spacing, starting from the con- 
tour of the sample edge. When there is a constriction, the 
stream lines must adapt by bending, which leads to for- 
mation of so called d-lines, i.e., lines in the flow pattern 
where the stream lines change direction discontinuously. 
These d-lines are recognized also in the corresponding 
magnetic flux distribution, due to an almost total local 
suppression of the magnetic field in their vicinity. By ex- 
amination of measured flux distributions, e.g., obtained 
by magneto-optical imaging (MOI), it has been found 
that the Bean model to a large extent explains the flux 
penetration patterns of hard superconductors. For ex- 
ample, the d-line created by a circular nonconducting 
hole (antidot) is parabolic [T5HT7] , a row of antidots give 
asymptotically straight d-lines at an angle given by the 
antidot fraction [14], while a convex right-angled corner, 
e.g., as in a rectangular sample, gives straight d-lines 
making 45° with the meeting edges [18] . 

Despite its vast success, the Bean model procedure 
for drawing critical current stream-line patterns is am- 
biguous at a concave corner. This is apparent in Fig. [I] 
demonstrating two possible ways to join stream lines as 
they pass a 90° concave corner. In (a) the equidistant 
stream lines are straight and meet at a 45° d-line, as in 
the case of a convex corner. In (b) the current stream 
lines consist partly of circular segments, and there are 
no traces of d-lines. The constructions (a) and (b) are 
just two out of infinitely many possible solutions that 
conserve current and give constant current density [19] • 
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(a) 




FIG. 1. Two different stream line patterns near a concave 
corner with flow of a constant-magnitude current density, j c , 
as in the Bean model. 



However, only one construction can be correct, and to 
identify it a more comprehensive analysis must be made. 

Schuster, Kuhn and Brandt [20] investigated the flux 
and current distributions in a cross-shaped superconduc- 
tor, a geometry that contains both convex and concave 
corners. Their MOI experiment showed no d-lines at the 
concave corners, hence excluding the sharp turn illus- 
trated in Fig. [T^i. An accompanying flux creep simula- 
tion gave stream lines close to circular, as in Fig. ^p. 
At the same time, a theoretical estimate showed that ex- 
actly circular stream lines diverge as E oc 1/r close to the 
corner and the i^-field will make a jump where the circu- 
lar stream lines connect with the straight lines running 
parallel to the edges. As pointed out by Ref. [21] this 
is an indication that the stream line pattern is unphys- 
ical. Thus, the stream lines at a concave corner cannot 
be perfectly circular. 

Gurevich and Friesen p~9j [21] calculated analytically 
the stream lines in the steady state using the powerful 
and general hodograph method. Their solution showed 
that, in the limit of negligible flux creep, neither Fig(l](a) 
nor (b) are correctly representing the stream line pattern 
near a concave corner. However, their solution did not in- 
clude the history of the magnetization process leading to 
the current distribution. Neither did it include the spe- 
cific non-locality effects crucially important in the trans- 
verse geometry. Thus, it still remains an open question 
what is the correct current stream line pattern in a thin 
superconductor with a concave corner, and subjected to 
a perpendicular magnetic field. 

This work considers in detail the flux penetration and 
current flow in a thin superconductor having the simplest 
shape that includes a concave corner, namely that shown 
in Fig. [I] A general simulation method for flux dynam- 
ics in thin samples of almost any shape is described in 
Sec. II. The approach takes special care to model the 
non-locality of the equations in an efficient way. Section 
III reports and discusses the simulations results which 
include distributions of magnetic field, shielding current 
and electrical field. A direct comparison with MOI ex- 
periments on a YBa2Cu30^ film in an increasing applied 



magnetic field is reported in Sec. IV. Finally, Sec. V gives 
the conclusions. 



II. SIMULATION MODEL 

The numerical scheme is quite general and can be ap- 
plied to thin planar superconductors of arbitrary shape, 
provided the thickness d is much smaller than any lateral 
dimension, and the whole sample can be embedded inside 
a rectangular area [14] . We also require that the external 
field H a is applied in the z-direction, transverse to the 
sample plane. Then, the flux dynamics is found by solv- 
ing Maxwell's equations in the quasi-static (eddy current) 
approximation, where the superconducting properties en- 
ters the equations through a highly nonlinear E — J rela- 
tions that characterizes the sharp vortex depinning tran- 
sition happening when the sheet current J approaches 
the value of the critical sheet current J c = dj c . A re- 
alistic approximation for many materials is a power law 

m 



E = pJ/d, p = Po ( J/J c ) r 



(1) 



where E is the electric field, p is the resistivity, and po 
is a resistivity constant. The creep exponent n is usually 
large, with the Bean model corresponding to the limit 
n — » oo. Note that the model also works for an ohmic 
conductor, where n = 1. 

Since current is conserved, V • J = 0, we can introduce 
the local magnetization g — g(r,t) as 



dg_ 
dy 



— Jxi 



dg _ _ T 

dx ~ v ' 



(2) 



where r = (x,y). Outside the sample, g = 0. The inte- 
gral of g gives the magnetic moment, m = J d 2 rg(r). 

For quasistatic situation, the B — J relation is given 
by the non-local Biot-Savart law, which can be rewritten 
to a B z — g relation as 



Bz/vo = H a + Qg, 



with the operator Q given by 
Qg{v) = J- 1 



where T is the 2D spatial Fourier transform, k 
k is the wave- vector. The inverse relation is 



<rV(r) = t- 



(3) 



(4) 



Ikl, and 



(5) 



where cp is an auxiliary function. Both Q and Q _1 are di- 
rect product in Fourier space and can thus be calculated 
effectively. 

By taking the time derivative of Eq. ([3| and inverting 
it, we get 
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B z /iiq — H Q 



(6) 



3 



This equation is solved by discrete integration forward in 
time. In order to carry out the time integration, B z must 
be known in the whole plane z = at time t and, in this 
work, completely different approaches will be used to find 
B z within the superconductor and in the surrounding 
vacuum. 

Starting with the superconductor itself, it obeys the 
material law, Eq. 0, which when combined with Fara- 
day's law, B z = — (V x E) 2 , gives 



B z = V.(pVg)/d. 



(?) 



From g(r,t) the gradient is readily calculated, and since 
the result allows finding J(r, t) from Eq. also p(r,i) 
is determined from Eq. ([I]) . The task then is to find B z in 
the non-superconducting parts, so that g = outside the 
superconductor. This cannot be calculated efficiently us- 
ing direct methods due to the nonlocal relation between 
B z and g. Instead we use an iterative procedure. 

For all iteration steps, i = B z ^ is fixed inside the 
superconductor by Eq. At i = 1, an initial guess is 
made for outside the sample, and #W i s calculated 
from Eq. (J6|. In general, this <jM does not vanish outside 
the superconductor, but an improvement can be obtained 

by 



(8) 



The projection operator O is unity outside the supercon- 
ductor and zero in the vacuum. Also, the output of the 
operation should be shifted to satisfy J cPrOg^ = 0. 
The constant is determined by requiring flux con- 
servation, 



d*r[B% 



(i+l) 



VoH a ] = 0. 



(9) 



Thus, at each iteration (i), B z l+1 ^ is calculated for the 
outside area. The procedure is repeated until after i = s 
iterations g^ becomes sufficiently uniform outside the 
sample. Then, g^ is inserted in Eq. (|6|, which brings us 
to the next time step, where the whole iterative procedure 
starts anew. 

The non-dimensional form of the equations are par- 
ticularly simple when the applied field is ramped with 
constant rate H a ^ and J c and n are both constants. 
We define the sheet current constant 



Jo = J c 



( dwp H a 
\ PqJc 



1/n 



(10) 



where w is some lateral length and the rest of the param- 
eters have been defined previously. The time constant is 
defined as 

t = J /H a . (11) 

The dimensionless form of the material law, Eq. 0, be- 
comes 



where J = J/ J and E = Ed/ (poJo)- The time evolution, 
Eq. (|6|), becomes 



dB z 
di 



(13) 



where g = g/(wJo), Q = Q/w, t = t/to, and B z = 
B z /{noJ(i). Finally, the Faraday law, Eq. becomes 



dB 



dt 



- = V- 



(14) 



where V = w\7. This means that the only free parameter 
in the model is n. A practical consequence of this is that 
simulations only need to be run once for each geometry 
and each value of n, and the result for any combination of 
dimensional parameters can be found simply by rescaling 
the solution. 

Several physical conclusions can be drawn directly 
from the non-dimensional equations. First, Eq.(13) with 



Eq. (14) inserted is a nonlocal diffusion equation with p 



E = pJ, 



jn— 1 



(12) 



as a very nonlinear diffusion constant. The non-locality 
prevents scaling solutions like in bulk [22 , but will still 
give plateaus in the current where J ~ 1. Second, the 
level of these plateaus in dimensional units is Jo, not J c 
as one could naively expect. Thus, J > 1 does not imply 
J > J c , since Jo < J c for parameters corresponding to 
most type-II superconductors. Third, the current con- 
stant depends on ramp rate as Jo ex H^ n '. Such a ramp 
rate dependent current gives also a flux penetration de- 
pending on ramp rate, and this has indeed been measured 
by Ref. [23] in a strip of high-T c superconductor. 

In this paper the tildes are omitted when reporting the 
results in dimensionless units. 

We chose the lateral length scale w as the half-width 
of the shortest side of the superconducting corner. The 
long sides of the sample have length 6. The sample is 
embedded in a 9 x 9 square, which is discretized on a 
512 x 512 grid. Note that since the geometry is very non- 
symmetric, only 25% of the grip points lay within the 
sample, the rest are in the vacuum, needed to fulfil the 
boundary conditions. Nevertheless, the method is fast 
and the simulations of this paper can easily be run on a 
personal computer. 



III. SIMULATION RESULTS 

The initial state is a flux free superconductor, and we 
use a flux creep exponent of n = 29. As the applied mag- 
netic field is increased with constant rate the magnetic 
flux gradually enters the sample from the edges, as seen 
in Fig. [2j showing snap-shots of the ^-distribution at 
H a = 0.3, 0.6 and 1.3. In the figure, the image bright- 
ness represents the magnitude of B Z1 which everywhere 
is directed parallel to the applied field. 

At H a = 0.3, the most visible feature is that the edges 
light up with field values much larger than the applied 




FIG. 2. Distribution of flux density B z at applied perpendic- 
ular fields H a — 0.3, 0.6 and 1.3. The images are color coded 
so that the brighter green the larger is B z , see the color bar. 
Included in each panel is the contour line of a constant flux 
density corresponding to the applied field, B z = noH a . 

field, while the sample interior is black in an area rep- 
resenting the Meissner state where the flux is expelled 
from the superconductor. Already at this low field, the 
concave corner shows enhanced field values and deeper 
flux penetration as compared to the penetration from 
the straight edges. 

At H a = 0.6, the flux front has advanced considerably, 
and now the concave corner is filled with a fan of en- 
hanced flux density. In the convex corners, the d-lines 
start developing, but they are not yet distinct lines, but 
rather wedge-shaped regions of vanishing flux density. 

At H a = 1.3, the flux has penetrated essentially the 
entire sample, and the superconductor is now described 
by a fully developed critical state. The d-lines are clearly 
seen in all convex corners, where they make 45° angles 
with the meeting edges, as expected for a superconduc- 
tor with isotropic J c . In the concave corner, on the other 
hand, there is no sign of any d-line. Quite the contrary, 
the flux density is there higher than in the surrounding 
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FIG. 3. Distribution of sheet current magnitude J and 
current stream lines at applied perpendicular fields H a = 0.3, 
0.6 and 1.3. The images are color coded so that the brighter 
red the larger is J, see the color bar. 



regions. In particular, in the corner point itself the mag- 
netic field is greatly amplified. Evidently, the current 
stream line pattern responsible for this flux distribution 
must be close to Fig.[lJ>, and definitely not like in Fig.[lji. 

Figure [3] shows color-coded images of the sheet current 
magnitude, J, together with the stream lines of J, cal- 
culated as the contour lines of the local magnetization g. 
At H a = 0.3 the critical state with J « 1 has only pene- 
trated a short distance from the edge, while most of the 
sample is in the Meissner state, where J < 1, as expected 
for a thin superconductor in a small magnetic field. The 
outermost stream lines can be followed around the sam- 
ple perimeter where they, along each straight edge seg- 
ment, stack with equal spacing, in accordance with the 
critical state model. At the concave corner the stream 
lines turn 90° in a gradual way, forming lines resembling 
circle segments while maintaining an essentially constant 
separation. Note that the two innermost stream lines, as 
they leave the critical-state in the concave corner, they 



FIG. 5. The distribution of electric field E close to the 
convex corner at H a = 1.3. 




FIG. 4. Distribution of electric field E at applied magnetic 
fields H a = 0.3, 0.6 and 1.3 during ramping of the magnetic 
field. The images are color coded to show the variation in 
the magnitude of E, see the color bar. Note that E is not 
calculated in the vacuum. 

cross over into the Meissner-state area. This is a clear 
manifestation of the non-locality of the governing equa- 
tions, and cannot be found within the bulk case of the 
Bean model. 

At H a = 0.6 the domains of constant current density 
have grown in size mainly by penetrating deeper, but 
also by filling larger parts of all corners. All the fea- 
tures commented for H a = 0.3 are here still present. In 
addition, the current now bends quite sharply at the con- 
vex corners, where the d-lines become increasingly more 
well-defined. 

At H a = 1.3, the full penetration state has been 
reached and the current density is everywhere J « 1. 
Note that in the central part of the sample, there are 
now several stream lines making small closed loops. At 
he same time, the d-line extending from the large con- 
vex corner is longer than the others. This overall pattern 
is fully compatible with the critical state model rule for 
constructing equidistant stream lines. In particular, one 



finds that the length of the long d-line should be a fac- 
tor 2(2 — a/2), or 17% longer than the others, in very 
good agreement with the simulation result. Moreover, 
the construction implies that the curved part of the cen- 
tral d-line consists of two parabolic curves meeting at the 
end of the long d-line. 

Figure [4] shows the distribution of the electric field 
magnitude E. At the straight edges, the results are as 
expected from the critical state model in an infinite strip, 
where the Enfield grows almost linearly from the flux 
front towards the edges [8]. In the Meissner state are 
E = 0. Also in agreement with the critical state model is 
the strong suppression of E close to the convex corners, 
where at the d-lines E = 0, signifying absence of any flux 
traffic. However, the most striking feature in the E-maps 
is the spot of very high field at the concave corner, which 
is strongly present in all three panels. Fig. [5] focuses at 
this spot at H a = 1.3. The results appears similar to the 
numerical result of Ref. [20], having a maximum value 
E ~ 4.5, compared to E ~ 1 at the straight edges. This 
enhanced electrical field is a sign of intensive flux traffic 
through the concave corner. 

Since high electrical fields are known to trigger ther- 
momagnetic avalanches [24-26], it is thus clear that thin 
superconductors with concave corners are far more sus- 
ceptible for such dramatic events to occur than samples 
with only convex corners or without corners at all. 

As mentioned, the current stream lines near the con- 
cave corner appears to be nearly circular. Yet, the sim- 
ulation results shown in Figs. [4] and [5] deviate from the 
E-fieid generated by perfectly circular stream lines, as 
calculated in Ref. [20 . The deviation manifests in three 
important ways. First, the field value at the edge is fi- 
nite, and is not showing the divergent E ~ 1/r behavior. 
Secondly, the maximum E is not found exactly at the 
corner, but in two spots located on each side of the cor- 
ner. Thus, the traffic of flux into the sample goes through 
a wider region than just a singular point, explaining why 
E is not divergent. Thirdly, the E-fie\d is smooth ev- 




FIG. 6. The current stream lines resulting from the simula- 
tion (full red line) and for comparison, stream lines including 
circular segments (dotted black line). 



erywhere inside the sample, and in particular there is 
no discontinuity in E at angles 0° and 90°. Hence, the 
stream lines obtained from the simulations are not sim- 
ply circular segments connected with straight lines. This 
motivates a closer inspection of the flow pattern. 

Shown in Fig. [6] is the detailed current stream line 
pattern near the concave corner at the applied field 
H a = 1.3. For direct comparison the figure shows also 
stream lines shaped as concentric circular segments con- 
necting to straight lines running parallel to the edges, and 
having a spacing corresponding to J = 1. The two pat- 
terns clearly deviate, as the calculated stream lines are 
compressed as they pass the corner, implying that there 
is an enhanced current density, J > 1, in the region. 
This is indeed consistent with the electrical field map in 
Fig.[5j where the highest electrical field is E ~ 4.5 giving 
a sheet current of 4.5 1 / 29 = 1.05, i.e., a 5% increase in 
the current density. 

Notice also that the calculated stream lines begin to 
curve some distance before they meet the corner sec- 
tor, thus the region of curves lines cover a sector wider 
than 90°. In this way, the stream lines change direc- 
tion more gradually, and the unphysical discontinuity 
in E mentioned earlier, is avoided. We also note that 
the stream line pattern in Fig. [6] deviates from the ana- 
lytically derived result presented in Ref. [19], where the 
stream lines make a sharp bend, near 45°, midway into 
the corner. Thus, it is evident that both the nonlinear- 
ity and nonlocality of the governing equations are impor- 
tant for the outcome of the modeling. One open question 
that remains is how the results would change in the limit 
n —> oo. The numerical calculation of this work does not 
give a definite answer, but since the applied creep expo- 
nent n = 29 is rather high, we expect that the results 
would be qualitatively the same. 



6 




FIG. 7. Magneto-optical images of flux penetration in a 
YBa2Cu3 x film at 45 K. The images in (a) and (b) were 
recorded at applied fields of B a = 16.5 mT and 44 mT, re- 
spectively. 



It is appropriate now to compare the simulations with 
results of MOI experiments. 



IV. EXPERIMENTS 

A film of YBa2Cu30 :c was made by laser ablation on 
a (100) SrTiOs substrate. Details of the preparation can 
be found in Ref. [27 . The sample has a thickness of 
d = 300 nm with the c-axis oriented perpendicular to 
the film plane. The critical temperature T c , measured by 
magnetic susceptibility, was 89.9 K. The critical current 
density of the film is j c = 8.0 x 10 10 A/m 2 at 45 K, 
the temperature where the magneto-optical images were 
recorded. 

The MOI investigation was performed using a bismuth 
substituted ferrite garnet film with in-plane magnetiza- 
tion as Faraday rotating sensor [28], placed directly on 
the sample surface. The sample was mounted on the cold 
finger of a continuous He flow cryostat with an optical 
window. Images of the flux distribution were recorded 
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with a digital camera through a polarized light micro- 
scope using crossed polarizers. In this way the image 
brightness represents the magnitude of the flux density. 
For more details of the method and setup, see Ref. [29] . 

Shown in Fig. [7] are two images of the flux penetration 
in the sample after an initial zero-field cooling to 45 K. 
Then the applied perpendicular magnetic field was slowly 
increased, and the images in (a) and (b) were recorded 
at B a of 16.5 and 44 mT, respectively. The image (a) 
was slightly contrast enhanced to allow the location of 
the flux front to be clearly seen. 

In the concave corner region it is evident that the flux 
front has advanced deeper than from the long straight 
edges of the sample. The swollen region covers a sector 
slightly exceeding 90°. The similarity between this ex- 
perimental image and the simulated result for H a = 0.6 
shown in Fig. [2] is striking. 

Also on the other side of the sample, near the main con- 
vex corner, the experimental images show a non-trivial 
behavior. The flux front associated with the lower edge 
penetrates slightly deeper near the main convex corner. 
The dashed line included in panel (b) is parallel to the 
edge, and serve as a guide to the eye. The formation of 
the extra long d-line in the main corner is also clearly 
seen in the image (b). The same swelling effect is seen 
to occur if one follows the similar flux front down along 
the vertical part of the sample. The agreement between 
these experimental results and the numerical simulations 



is striking down to the very fine details. 

V. CONCLUSIONS 

In superconducting films, the electrodynamic non- 
locality implies that calculations, in principle, must con- 
sider the whole sample and the full magnetic history. 
This is particularly evident at the concave corner consid- 
ered in this work, where the current stream line pattern is 
nontrivial, unlike the convex corner where the stream line 
pattern can be drawn using a simple Bean model proce- 
dure. The flux creep simulation conducted in this work 
shows that the stream line pattern at the concave cor- 
ner deviates from exact circles, and this small deviation 
prevents unphysical jumps and infinities in the electrical 
field. Nevertheless, the electrical field is very high at the 
concave corner, and samples with such corners should 
thus be particularly susceptible for nucleation of thermo- 
magnetic avalanches. The magneto-optical experiment 
of this work shows striking similarity with the simula- 
tion and verifies the necessity of properly including the 
electrodynamic non-locality when modeling thin films in 
transverse geometry. 
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